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LARGE TIME-STEP STABILITY OF EXPLICIT 
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Lewis Research Center 
Cleveland, Ohio 44135 

and The University of Akron 
Akron, Ohio 44325 

ABSTRACT 

There is a wide-spread belief that most explicit one-dimensional advection schemes need to 
satisfy the so-called “CFL condition” — that the Courant number, c = uAt/Ax, must be less 
than or equal to one, for stability in the von Neumann sense. This puts severe limitations 
on the time-step in high-speed, fine-grid calculations and is an impetus for the development 
of implicit schemes, which often require less restrictive time-step conditions for stability , 
but are more expensive per time-step. However, it turns out that, at least in one dimension, 
if explicit schemes are formulated in a consistent flux-based conservative finite-volume form, 
von Neumann stability analysis does not place any restriction on the allowable Courant 
number. Any explicit scheme that is stable for c < 1, with a complex amplitude ratio, G(c), 
can be easily extended to arbitrarily large c. The complex amplitude ratio is then given by 
exp (~iN6) G(Ac), where N is the integer part of c, and Ac = c - N (< 1); this is clearly 
stable. The CFL condition is, in fact, not a stability condition at all, but, rather, a “range 
restriction” on the “pieces” in a piece- wise polynomial interpolation. When a global view 
is taken of the interpolation, the need for a CFL condition evaporates. A number of well- 
known explicit advection schemes are considered and thus extended to large At. The 
analysis also includes a simple interpretation of (large- Ar) TVD constraints. 

INTRODUCTION 

Consider the model one-dimensional pure advection equation for a scalar <t>(x,t) 

<>± = _ u (i) 

dt dx 

where u is a constant advecting velocity. Take a uniform space-time grid (Ax, Ar) and 
integrate Equation (1) over Ax and At, giving a finite-volume formulation 

if = i; - c «>, - *,) (2) 

where the bars indicate spatial averages over cell i at time-levels n and (n+1), and left and 
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right time-averaged face-values have been introduced. The Courant number is given by 

c = uAt/Ax (3) 

The notation is defined in Figure 1, which also shows advective characteristics (along any 
one of which 4> is constant) entering the left face (for u > 0). Given the set of current 
spatial-average values, <£", one needs to estimate the corresponding time-averaged face- 
values in order to explicitly update according to Equation (2). 

In a von Neumann stability analysis, <f>(x,t) is written as a wave [Fletcher, 1990] 

4>(x,t) = -4(0 exp(t kx) < 4 ) 

where k is the wave number and t represents the imaginary unit, ■ When this is 
substituted into Equation (1), the exact solution is [Leonard, 1980] 

<t>(x,t) = A(0) exp[i k(x-ut)] @) 

corresponding to a travelling wave, with <t> = const along characteristics (x = ut). The 
exact complex amplitude ratio is then 

G = <£ (*_, i+AQ = (6) 

e “ ct 4>(x,t) 

where 6 is the nondimensional wave-number 

6 = kAx V) 

Note that if c > 1 , G exiCt can be written 


where N is the integer part of c 


= exp(-t N6) exp(-t Ac 6) 


N = INT(c) 


and Ac the non-integer part (necessarily less than one) 

Ac = c-N < 10) 

By integrating Equation (4) from (x, - Ax/2) to (x ; + Ax/2), it is not difficult to show 
that the exact cell-average values are 


(/>" = .4(0) exp{i ^[^ -w/zAt]} 


sin 8/2 


</>" M = A(0) expJifct^.-uOi+OAf]} 


sin 61 2 
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Hence, the exact complex amplitude ratio of cell-average values is the same as that given by 
Equation (8) 

T»+l 

G cxict = _j— = G ejuct = exp(-tlW) exp(-iAc0) ( 13 ) 

thus, giving a reference value against which to compare numerical G’s. 

In the following sections, the von Neumann stability of some well-known explicit 
advection schemes is studied — first in terms of the “conventional” time-step ranges, and 
then for larger Courant number values. The latter extension stems naturally from identifying 
the sub-grid interpolation used in estimating the time-averaged face-values in Equation (2). 
It will be shown, in particular, that if Gf^(c) is the (conventional, small time-step) complex 
amplitude ratio for a given scheme for c ^ 1, then the natural finite-volume extension to 
arbitrarily large time-step has a complex amplitude ratio of the form 

i 

G ™(c) = ZL- = exp(-.W) G™(Ac) for Ol (14) 

which should be compared with Equation (13). To the extent that q*™(Ac) is a good 
approximation to exp(-iAc0), then the large- Af G^(<0 is an even better approximation 
to the exact G. Viewed differently, a numerical simulation (over a given distance) will be 
more accurate for larger A t because, as will become evident, numerical distortion depends 
only on the total number of time steps and Ac — not on c itself. 

EXPLICIT ADVECTION SCHEMES 


The time-averaged face- values appearing in Equation (2) can be rewritten in terms of spatial 
averages. For example, the time averaged left-face value is given by 


<t>i ~ 



1 

c Ax 


<*>"(£) 


' x, -cAx 


(15) 


where ( T ) is the instantaneous face-value and </>"(£) is the local sub-grid spatial behaviour 
in the region upstream of the face in question, at time-level n. The relationship should be 
clear from Figure 1; note that uAt = cAx. A similar formula holds for the right face. But 
it is not necessary to write this out explicitly because advective flux conservation guarantees 


that 


4> r (j) = <M I+1 ) 


(16) 


Different numerical schemes result from different choices of <£"(£) in estimating the local 
behaviour. Note in particular that for consistency, <£"(£) should obey the integral constraint 


— = 4>; for all i 

Ax l xr ^n 


( 17 ) 
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First-Order Vpwinding 

One of the simplest advection schemes results from assuming that <£"(£) is piece- wise 
constant, for each i, 

*"({) = hi for (*,-y) <5 < (^y) <18> 

with discontinuities at the cell faces, as shown in Figure 2, the hatched region depicting the 
integral in Equation (15). This, of course, trivially satisfies Equation (17), and the left-face 
value given by Equation (15) for a positive Courant number less than (or equal to) one is 
simply 


4> t (i) = 4>U for 

0 < c < 1 

(19) 

and, similarly 

<t> r ( i) = 4>* for 

0 <c < 1 

(20) 

The update equation, Equation (2), thus becomes 

Jr 1 = 5; - c ti] - *U) 

for 0 < c < 1 

(21) 

The corresponding complex amplitude ratio is 

G w (c) = 1 - c[l 

- exp(-*0)] 

(22) 

or 

Giu(c) = 1 - c(l-cos0) - icsinfl 

(23) 


For the full range of numerical wave-numbers (0 < 6 < t), this represents a semicircle of 
radius c in the lower half of the complex plane, passing through the point (1,0). The scheme 
is thus stable (|G| < 1 ) provided the CFL condition ( c < 1) is satisfied. Figure 3 shows 
a polar plot of G m for c = 0.75. 

Making a Taylor expansion of real and imaginary parts gives 

G 10 (c) = 1 - £ 0 2 + 0(6*) -i (cO + O(0 3 )] (24) 

whereas a Taylor expansion of the exact G gives 

G exact (c) = cos(c0) - t sin(c0) 

= 1 - Sle 2 * 0(6*) - t [cd - Si 6 3 + O(0 5 )] (25) 

Thus Gju indeed matches G ejuct only through first-order terms in 6, with discrepancies in the 
second-order term (unless c = 1 , in which case exact cell-to-cell transfer occurs across one 
mesh-width: ^"* 1 = <^" 1 ). 


i 
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Second-Order Methods 


It is convenient to define the “order” of an advection scheme as the power of 0 up to and 
including which the Taylor expansion of G num matches that of G exict given by Equation (25). 
One of the best known second-order methods is that due to Lax and Wendroff [1960], 
equivalent in the scalar case to Leith’s method [1965]. For this scheme, the sub-grid <£"(£) 
is assumed to be piece-wise linear satisfying Equation ^17). Upstream of the left face, for 
example, a straight line is drawn (for u > 0) between and <f >" » considered to be located 
at the centers of the respective finite-volume cells. A similar construction is made across 
each of the other cells, resulting in discontinuities at cell faces. This is shown in Figure 4, 
which also indicates the integral in Equation (15) for 0 < c < 1. 

In this case, i.e., for a Courant number less than (or equal to) one, 


<*>"(£) 



for 




<£ < 



(26) 


upstream of the left face, across cell (i-l). Substitution into Equation (15) results in 


<M0 




(27) 


According to Equation (16), <t> r is obtained by increasing each index by 1. This gives the 
update equation as 






- 2 <}> i + ^..j) 


(28) 


Because of the symmetry about cell i, this formula is actually valid for positive or negative 
Courant number, provided 

|c| * 1 < 29 > 

The negative-c case is shown in Figure 5, where the sub-grid 4> n (£) for computing 4> t is the 
same formula as that in Equation (26) — but across cell i rather than (t-1). 

The numerical complex amplitude ratio for the Lax- Wendroff scheme can be obtained 
directly from Equation (28) as 

G lw (c) = 1 - £[exp(i 0) - exp(-i0)] + y[exp(i0) - 2 + exp(-i0)] ( 30 ) 


or 


G (c) = 1 - c 2 (l - cos0) - ic sin 0 


with a corresponding Taylor expansion 

= 1 - Si 6 2 + <7(0 4 ) - t [c6 - £ 0 3 + <9(0 5 )] 


(31) 


G lw (c) 


(32) 
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thus matching the expansion of G exact , Equation (25), through 0 2 terms, as expected for a 
second-order scheme. Once again, exact point-to-point transfer across one mesh-width 
occurs when c = 1 (and also, in this case, when c = -1). Figure 6 shows a polar plot of 
G lw , a semi-ellipse with vertical semi-axis equal to c and curvature near (1,0) equal to that 
of the unit circle. The scheme is stable provided the CFL condition is satisfied — and 
unstable otherwise. 

The Lax-Wendroff scheme is symmetrical in the sense that a given face-value 
(assuming | c | < 1) depends only on the two immediately adjacent cell-average values on 
either side of the face in question. Figure 7 shows an alternate piece- wise linear sub-grid 
reconstruction (for positive c) using the two immediate upstream cell-average values. 
Although the reconstruction of 4>*(£) looks superficially the same as that in Figure 5, this 
represents an entirely different advection scheme in that the integral of Equation (15) is 
computed from the left (rather than the right, as in Figure 5). This will be found to generate 
a second-order advection scheme. Because of the upwind bias, it is now commonly known 
as “second-order up winding” [Fletcher, 1990]. Upstream of the left face, across cell (i -1), 
</>"(£) is given by 


<*>"(£) = + 


Ax 


(£ for 


Equation (15) then gives the left-face value as 

♦. - M *» - tg) 


(*■-• - 


Ax 

2 


< I < 


K* 


for 0 < c < 1 


Ax | (33) 


(34) 


In fact, this equation is valid for 0 < c < 2. This is an important point (usually ignored 
in CFD literature) because, as described below, it suggests a natural and computationally 
efficient way of extending explicit advection schemes to arbitrarily large Courant numbers. 
Using the corresponding 4 > r , the update equation for second-order upwinding becomes 


4>t = 4>t ~ c 


- (2 - c) K, 

and the complex amplitude ratio is 

G 2U (c) = 1 - c (2^£) + c (2 - c) cos0 - c cos 20 


- t c(2 -c) sin0 - c( LJL ) sin20 


(35) 


(36) 


This has a Taylor expansion given by the following 
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g 2U (c) = i - c - 0 2 + o<0 4 ) - * 


c0 - ( 3c2 — ) 0 3 + 0(0 5 ) 


(37) 


Note in particular that, for c = 1, 


G 2U (1) = cos 0 -isin0 = G e3UCt (l) (38) 

corresponding to exact cell-to-cell transfer across one mesh- width. But, in addition, note 
that for c = 2, 

G 2u ( 2) = cos 20 - i sin 20 = G exact (2) (39) 

with exact cell-to-cell transfer across two mesh-widths. In fact, second-order upwinding is 
stable throughout the entire range 


|G 2U | <1 for 0 < c < 2 ( 4 °) 

Two typical polar plots (for c = 0.75 and 1.75) are shown in Figures 8 and 9, respectively. 

In order to prove (40), it is convenient to rewrite G 2U (c) in the following form, 
starting with Equation (36), 

G 2u (c) = (cos0 - i sin 0) [ 1 - (c - l) 2 (1 - cos 0) - i (c - 1) sin0 ] (41) 

Note, by comparing Equation (31), that this is of the form 


G 2U (c) = exp(-i0) G lw (Ac) 


(42) 


where Ac = c - 1. Of course, |exp( -i0) | = 1, and for the “Lax-Wendroff” component 
based on Ac, 

| G lw (Ac) | <1 for -1 < Ac < 1 ( 4 3) 


thus proving (40). It should be stressed that second-order upwinding is one well-known 
explicit advection scheme for which the CFL condition does not apply. 

In order to gain some insight into the operation of second-order upwinding for a 
Courant number range of 1 < c < 2, Figure 10 shows the situation for c = 1.25. In this 
case, it is not hard to show, using Equation (15), that 


c<t> t 


Ax 


x ( -&x 

{ <t> n (Z)dZ 

x,-c Ax 


JL_ 

Ax 


f 

J x,-*x 


= ac (4>u * *;. 2 > - ^ - *;.o] + c. w 

with a similar formula for c<t> r (with all indexes increased by 1). The term in square 
brackets is clearly the Lax-Wendroff face- value based on Ac — but for cell (i-1) rather than 
cell i. This means that the update equation takes the form 



8 



<*>;., - ^ «; - 


Ac 2 


(<}>: - 20"_j + *;_ 2 ) 


(45) 


with Ac = c - 1. This is algebraically identical to Equation (35). But, by comparison with 
Equation (28), it can be seen that the right-hand side is the Lax-Wendroff update of cell 
(i-1) using Ac rather than c itself. This gives a hint as to how flux-based conservative 
finite-volume explicit advection schemes might be extended to arbitrarily large At. The 
details will be examined in the next section. 

As is well known, the conventional Lax-Wendroff scheme suffers from severe phase- 
lag dispersion [Fletcher, 1990]. On the other hand, second-order upwinding introduces 
phase-lead dispersion for 0 < c < 1. Fromm’s method [1968] of “zero-average-phase- 
error” attempts to minimize the dispersion by taking the arithmetic-mean of the respective 
sub-grid reconstructions, giving, from Equations (26) and (33) 


4> B (£) 



,n ,n 

-4>i- 2 \ 
l 2 Ax I 


(l for (*i- 1 


Ax 

~2 


) <£< 



(46) 


i.e., the slope across any cell is parallel to the chord joining <£ -values through the centers 
of the two adjoining cells on either side. This is shown in Figure 11. To compute the face- 
value, first rewrite Equation (34) for second-order upwinding as 


4> t (2U) = I (</>; + 4>U) - \ to - O - W - 2 <t>U + 4>u) 


(47) 


which is seen to be the Lax-Wendroff value together with a correction term proportional to 
the upwind-biased second-difference. Then the left face-value for Fromm’s method is seen 
to be, for Courant number less than (or equal to) one, 

4>, (Fromm) = I [tf>,(LW) + <£,( 2U)] 


2 v ‘ ' 

and the corresponding update equation becomes 


= \ to * \ w - *u) - (L^) to - 2*;., * o 






with a complex amplitude ratio 


G Fr (c) = G lw (c) + c(-L-£) [exp(i0) - 3 + 3exp(-i0) - exp(-i20)] 


(48) 


(49) 


(50) 


or 
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G Fr (c) = 1 - c 2 (l - cos 0) - c(i— -) (1 - cos 0) 2 

(51) 

- i [c sin 0 + c{ l ~ ) sin 0(1- cos 0)] 

This has a Taylor expansion 

G Pr (c) = 1 - C -1 P * 0(6') - . [c# - * 0(9 ! )] (52) 

Although this is only formally second-order accurate, the 0 3 coefficient is a fairly good 
approximation of the exact coefficient, c 3 / 6, over the operating range, 0 < c < 1. From 
its formulation as the average of Lax-Wendroff and second-order upwinding, Fromm s 
method clearly needs to satisfy the CFL condition for stability. 


Higher-Order Methods 


In order to construct a higher-order approximation for estimating a face- value, say <j > t , one 
might try to interpolate a parabola across cell (i-1) passing through the three values, 
4> n , and 4> n (for u > 0), with the location of these values considered to be at the 
center of their respective cells. Such a parabola is shown by a dashed line in Figure 12. 
But, of course, this does not generally satisfy the integral condition given by Equation (17). 
Because of its curvature, the parabola needs to be shifted vertically to satisfy this constraint. 
Thus, by introducing an additional constant, <£"(£) can be written, in the upstream vicinity 
of the left face, as 


<*>"(£) = C, + 4>l 

fOT Z S (*M * 

then gives 








(Z-xJ 


(53) 


Ax\ 


Tax /" \ 2 Ax 2 

Substitution into Equation (17), rewritten for cell (i-1), 


C x = - 


24 


(54) 


The corresponding parabolic segment across cell (i-1) is shown by a heavy curve in Figure 
12. Then Equation (15) gives the corresponding left-face value, for Courant number less 
than (or equal to) one, as 


*, = \ (*; + - § (*: - +;-) - (i^) (*: - + *«) 


(55) 


with a similar formula for <f> r (i - i+1, as usual). The corresponding update equation is 
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ir = 1 (c 2 *;♦*•,) 

* c (*<•■ “ 3 *i * 3 *"-i ‘ *"-j I (56) 


This has the same form as the third-order upwind QUICKEST scheme [Leonard 1979], but 
in terms of cell-average values, ^ , rather than nodal values. But note from Equations (53) 
and (54) that a nodal value can be retrieved from the sub-grid behaviour. For cell (i-1), for 
example, the nodal value is 


*„ = *■(*,.) = 5;., - 


(57) 


This means that the QUICKEST scheme can be constructed from Equation (56) by writing 




*♦1 ,«♦ 1 / 
1 = *> "l 


, * + l ,l»+l 

<Pi+ 1 ~ 2< Pi +< t > i - 1 

24 


for all i 


(58) 


giving 


<t>r = j (Ci -*;) + y (c, -2ti+4>u} 

+ c (^T~) K 1 " + " ti- 2 ) 


(59) 


for 0 < c < 1, thus matching the original QUICKEST algorithm, for pure convection. 

From the similarity with Fromm’s method, the complex amplitude ratio can be 
immediately written down as 

G q (c) = 1 - c 2 (1 -cos 6) - c | 1 c ) (1 -cos 6 f 


-1 c sin0 + c(i — sin0 (l-cos0)j 


(60) 


and, in this case, the Taylor expansion matches G eMCt , Equation (25), through third-order 
terms 

G q (c) = 1 - — e 2 + O( 0 4 ) - t [c6 - sl 0 3 + 0(6 5 )] (61) 


verifying the expected third-order accuracy. 

Higher-order methods can be constructed in a similar fashion. If a Pth-order piece- 
wise polynomial (or spline, for that matter) is used for the sub-grid reconstruction, satisfying 
Equation (17), the resulting advection scheme will be (P+ l)th-order accurate. This pattern 
has already been seen for first- through third-order methods. 


it 
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LARGE TIME-STEP EXTENSION 


Most conventional explicit advection schemes are constructed on the assumption that the 
Courant number will be in the range 0 < c < 1 in computing time-averaged face-values; 
note the conditional statements just before Equations (19), (26), (48), and (55). It is not 
surprising, then, that such methods need to satisfy the CFL condition (|c| < 1) for stability 
— the notable exception being second-order upwinding, which, somewhat fortuitously, is 
stable for 0 < c < 2. Referring again to Figure 10, it can be seen that second-order 
upwinding, for c > 1, is equivalent to a “Lax-Wendroff type” of sub-grid reconstruction 
over each cell (compare cell i - 2 in Figures 4 and 10); i.e., for each cell i, 


*"(€) 



(62) 


for u > 0. Then, since the integration in Equation (15) is over more than one cell, the 
average left face, for example, can be written 


c*,(i) = 4>U + Ac <^ w (i-l, Ac) ( 63 ) 

where the notation implies that <^ w (i-l, Ac) is the left-face value of cell (z'-l) that would 
have been obtained for a Lax-Wendroff scheme with a Courant number of Ac (= c - 1, in 
this case). 

Clearly, this idea could be extended to even larger Courant number. For example, 
Figure 13 shows a situation for c = 2.625. In this case, it should be clear that 


c<t> t (i) = 4>U * 4>U + Ac^ w (i~2,Ac) 
where Ac = c - 2. Of course, the right-face value is given by Equation (16), or 


(64) 


c<t> r (i) = 4>" + + Actf| W (i-l, Ac) (65) 

The update equation then becomes 

fi’" = -[£♦;*,♦ Ac # w (/-l, Ac)] Ac*r(>-2,M (66) 

resulting in an effective update of the form 


= iu - Ac K W (i-l. Ac) - ti’V- 2, Ac)] (67) 

or, more specifically, using the sub-grid (i.e., “Lax-Wendroff-type”) reconstruction of 
Equation (62), 




Ac 2 

IT 


(<t>U -UU + 4>u) 


( 68 ) 


It is not hard to see that this idea generalizes to (in principle) arbitrarily large Courant 
number, in which case 
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c4> t {i) = ti-p + Ac Ac) 

p-1 

with a corresponding formula for c<f> r (i), according to Equation (16). 
(“extended Lax-Wendroff”) update equation is then 


(69) 

The effective 




Ac n 

2 


(^i-jV+l ~ ^i-N-l) + ~2~ fai-N* 1 ^ + ^i-N- 1) 


(70) 


Of course, the same principle could be applied to any sub-grid reconstruction. Then, 
one simply writes for the left face, say, 


c4> t (i) 


E *U 


+ Ac <f> t (i-N, Ac) 


(71) 


where 4> t (i-N, Ac) is computed at the left-face of cell ( i-N) based on a Courant number of 
Ac. The update equation is then 


<t> T 1 = </>;- [c^> f (i + l) - c<t> t (i)] 


(72) 


as usual. And this is equivalent to 

~4>r = * { . N - Ac W'ii-N+hAc) - <f> t (i-N,Ac)] 


(73) 


Stability Analysis 

Consider a von Neumann stability analysis of the extended Lax-Wendroff update equation, 
Equation (70). The complex amplitude ratio is 


G elw = exp (-iN0) - ^ jexp[-i(N-l)0] - exp[-i(N+l)0]} 


A c : 


{exp[-i(N-l)0] - 2 exp(-tM) + exp[-t(AT+l)0]} ( 74 ) 


which can be rewritten as 

G elw = exp(-iiW) {l - ^ [exp(t 6 ) - exp(-»0)l 

+ [exp(i0) - 2 + exp(-i0)]j 

or, by comparison with Equation (30), 

G elw = exp(-t NO) G lw (Ac) 


(75) 


(76) 
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where Ac = c-N. Since Ac < 1, by definition, the extended Lax-Wendroff scheme is 
clearly stable. 

The formula can be generalized to any scheme (i.e., any sub-grid reconstruction) for 
which the face- value is a linear combination of cell-average values. For such schemes, the 
large time-step (LTS) G can always be factored into the form 

G lts = exp (-.«) G cts (Ac) (77) 

And if |G sts (Ac)| < 1, for Ac < 1, the LTS scheme is stable. The formal order of 
accuracy of G LTS is the same as that of G STS (Ac). 


LARGE TIME-STEP EXPLICIT TVD SCHEMES 


Total-variation-diminishing (TVD) schemes, as described by Sweby [1984], for example, 
involve face-values that are, in general, nonlinear functions of (local) cell-average values. 
The TVD constraint boundaries described by Sweby have a simple interpretation in terms 
of sub-grid reconstruction. The explicit large- time- step extension is then quite straight- 
forward, although the von Neumann stability analysis is somewhat more involved (even for 
the case of c < 1). 

First, for definiteness, consider the case of c < 1. The left-face value is based on 
a linear reconstruction across cell (z’-l) satisfying Equation (17), but with a slope determined 
by certain “limiter” constraints. It is convenient to define 


<t>, = 0 ) 

This is shown in Figure 14. The TVD value of is usually written as 


<*>?( TVD) = ~ <P [*m - \ W + 0] 


(78) 


(79) 


interpreted as first-order upwinding, corrected by a term proportional to the difference 
between first-order upwinding and second-order central (i.e., linear) interpolation between 
4> n j and 4>" > as shown in the figure. The factor <p is the antidiffusion flux-limiter factor, 
and is a nonlinear function of the local gradient ratio 


r 


(♦‘-I ~^m) 


Equation (79) can be rewritten as 


<£°(TVD) 




(80) 


( 81 ) 


The time-averaged face-value (for c < 1) is then, from Equation (15), 
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giving 


TVD) = (1 -c)4>] + c<t>U 


<t>'( TVD) = <^ w - 1 (1 -*>)(! ~c) 


where the coefficient of the latter term can be viewed as being proportional to an effective 
artificial diffusivity (or viscosity) 

i (1 -<p)(l - c)kAx < 84 > 

relative to the Lax-Wendroff face-value. 

Sweby gives the TVD constraints on <p as 

<p = 0 for r < 0 ^ 8 ^ 

0 < <p < 2r for 0 < r < 1 ( 86) 


0 < <p < 2 for r > 2 

This is shown in Figure 15. Figure 16 shows three typical examples for r = - 1/2, + 1/3, 
and 3. For a local maximum (or minimum), case (a), 

<t>< = <t>°t = iU , r < 0 (88) 

In the other cases, the sub-grid reconstruction across cell (i - 1) is a linear function passing 
through at x, 

*’«) = ♦*«-*<-•) (89) 

where the slope, 5, is constrained to satisfy 

0 < |S| S |S„| (90) 

For a gradient increasing (in magnitude) with i (i.e., 0 < r < 1), case (b), the constrained 
slope is 

_ (ih-iU) (91) 

con Ax/2 

corresponding to <p = 2r, or = 2^ - <£" 2 . And for r > 1, case (c), 

C = 

con Ax/2 


( 92 ) 


15 


corresponding to <p = 2, or = <f>" ■ The limiting slope, S con , is thus seen to be the 
smaller (in magnitude) of twice the backward or forward gradient; the larger value is shown 
by the chain-dashed line in each case in Figure 16. 

Although explicit TVD schemes have been used only for time-steps satisfying the CFL 
condition, there is clearly no need for this constraint. Figure 17 shows a reconstruction that 
uses S con across each cell, with a Courant number of c — 2.25. In this case, 


c 4> t (i) = </>"_! + 0". 2 + Ac <f> t (i-2) 


(93) 


where 

Ac <*>,(/ -2) = Ac[#_ 3 + (1 - Ac) S coo Ax/2] 
or, in this case (since r > 1), 

Ac<£,(j-2) = Ac ^"_ 3 + Ac (1 - Ac) (4>1 2 ~4>U) 

And the right-face value is given by 

c<£ r (/) - c<f> ( (i+ 1) = <£" + tf-i + Ac<£,(z-1) 
or (since, here, r < 1) 

c 4> r (i) = # + i-.j + Ac 4>*. 2 + Ac (1 - Ac) (S"_ 2 - 4>U) 
for the case shown. The update equation is then, effectively, 


(94) 

(95) 

(96) 

(97) 


in this case, equivalent to 
(for cell 0 is 

G, 


*r = 


(98) 


(extended) first-order up winding. The complex amplitude ratio 
= exp(-i20) {1-0.25 [l-exp(-*0)]} (") 


which is stable, the term in square brackets being G 10 (Ac). 

In general, for c = N + Ac, the complex amplitude ratio will vary from cell to cell, 
depending on the particular form of sub-grid reconstruction near cell ( i-N ), but, in any 
case, 

G, = exp(-iM) G._ JV (Ac) ( 10 °) 

using an obvious notation. Since, for TVD schemes, | G^ N (Ac) \ < 1 f° r Ac < 1, the large 
time-step extension is always stable. 
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CONCLUSION 


Most conventional explicit advection schemes need to satisfy the CFL condition (c < 1) 
because of a restriction on the range of validity of the sub-grid reconstruction formula 
upstream of individual finite-volume faces. If the given formula for the sub-grid 
reconstruction (valid for c < 1) is extended beyond c = 1, the scheme is unstable in the von 
Neumann sense (except for second-order upwinding, which is stable for 0 < c < 2). For 
example. Figure 18 shows an extension of the conventional Lax-Wendroff sub-grid 
reconstruction upstream of the left face of cell i , for c = 1.75, 


*"(?) 




( 101 ) 


Clearly, the real problem is that, across cell (/ - 2), <£"(0 does not satisfy the consistency 
relationship of Equation (17); i.e., 


f r(t)dz * i“ 2 < 102 ) 

X,i-Axl2 

In the case shown, the influx through the left face of cell i will be too jmall (and the outflux 
through the right, far too large) resulting in a wild oscillation of <f>. that rapidly grows 
without limit. 

Second-order upwind sub-grid reconstruction upstream of the left face of cell /, 
(somewhat fortuitously) satisfies Equation (17) for both cell (i - 1) and (i -2). This is the 
fundamental reason for the stability range extending to c = 2. Once this is recognized, it is 
not hard to see how to generalize to arbitrarily large Courant number. One simply performs 
a sub-grid reconstruction satisfying Equation (17), then computes all face- values according 
to Equation (15), or, equivalently, 


C 4> t (i) = Y, 4 1-p + Ac <j> ( (i-N,Ac) 

p‘l 

(103) 

finally updating via 



(104) 

The result is equivalent to 



C = - Ac[*,(t-N + l,Ac) - <f> t (i-N,Ac)] ( 105 ) 

representing an update of cell-average 2 for a Courant number of Ac, together with a 
simple translation over N mesh- widths. It should be clear that any numerical distortion of 
the evolved profile depends on the particular form of sub-grid reconstruction, the total 
number of time-steps, and Ac— but not on c itself. 

A piece-wise constant assumption across each cell results in a very diffusive first-order 
scheme. A Lax-Wendroff-type sub-grid reconstruction, as given by Equation (26) for each 
cell, gives rise to trailing oscillations, just as in the conventional (c < 1) case. Conversely, 


n 
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the large time-step (LTS) extension of second-order upwinding, given by Equation (33) for 
each cell, results in leading oscillations. The extended version of Fromm s method, 
Equation (46), and particularly QUICKEST, Equation (53), are much less oscillatory. 
Finally, for large- A? explicit TVD schemes, the evolved shape is nonoscillatory, identical 
to that for a Courant number of Ac (< 1) for the same number of time-steps, but translated 
N mesh- widths downstream. 

Note that some sub-grid interpolations are independent of the direction of the 
convecting velocity. This is the case, for example, with (the large-At extension of) first- 
order upwinding, Fromm’s method, and QUICKEST’ spiecewise-parabolic interpolation; the 
reconstruction is the same, no matter whether c is positive or negative. In these cases, there 
is a sort of built-in “natural upwind bias’’ because of the way the fluxes are calculated. A 
second-order upwind type reconstruction based on (for c > 0) 

= ♦; * (*#) « - *,> (106) 

across each cell i, has an additional upwind bias. For c < 0, the corresponding formula 
would be 


4>\Z) 



(107) 


By contrast, the Lax-Wendroff type reconstruction uses Equation (107) for c > 0 and 
Equation (106) for c < 0; this is a form of downwind bias. Generally, the direction- 
insensitive schemes have better phase accuracy than the “directional” schemes. 

A von Neumann stability analysis of the large-time-step explicit scheme always results 
in a complex amplitude ratio of the form 

G lts = exp (-,N6) G^A c) (108) 


which should be compared with the exact expression, rewritten here for convenience 

G exic , = exp ( - iN6) exp ( - 1 Ac 0) (109) 


Generalization of the one-dimensional LTS schemes to variable mesh-widths and 
(some specific forms of) varying velocity fields appears to be relatively straight-forward. 
Such schemes have recently been explored by Roache [1992]. Extension to two (and three) 
dimensions is much more challenging; but some progress has been made. This will be 
reported in future publications. 
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Figure 11 . — Sub-grid reconstruction for Fromm's method. The slope of the interpolant across 
any cell is parallel to the chord joining adjacent cell -averages. In this case, c = 0.75. 



Figure 12. — Piece-wise parabolic interpolation across cells, 
scheme equivalent to QUICKEST; c = 0.75. 
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